1 Packages

library(ggplot2) #for plotting
library(tidyverse)
## ── Attaching packages ────────────────────────────────────────────────────────────────────────────────────────── tidyverse 1.3.0 ──
## ✓ tibble  3.0.1     ✓ dplyr   1.0.0
## ✓ tidyr   1.0.2     ✓ stringr 1.4.0
## ✓ readr   1.3.1     ✓ forcats 0.5.0
## ✓ purrr   0.3.4
## ── Conflicts ───────────────────────────────────────────────────────────────────────────────────────────── tidyverse_conflicts() ──
## x dplyr::filter() masks stats::filter()
## x dplyr::lag()    masks stats::lag()
library(dplyr)
library(GGally)
## Registered S3 method overwritten by 'GGally':
##   method from   
##   +.gg   ggplot2
library(gridExtra)
## 
## Attaching package: 'gridExtra'
## The following object is masked from 'package:dplyr':
## 
##     combine

2 Read in Data

The Pinery dataset was compiled in excel from monthly mean temperature and total precipitation data using open access Environment Canada records for stations within 30km of the Pinery with at least 10 consecutive years of data (http://climate.weather.gc.ca/historical_data/search_historic_data_e.html).

The London datasets are homogenized data from Environment Canada that were converted to tidy format in excel.

pinery <- read.csv("data/pinery_climate.csv")
pinery$station <- as.factor(pinery$station) #Make a factor 
str(pinery)
## 'data.frame':    2931 obs. of  5 variables:
##  $ year        : int  1882 1882 1882 1882 1882 1882 1883 1883 1883 1883 ...
##  $ month       : int  7 8 9 10 11 12 4 5 6 7 ...
##  $ mean_temp   : num  NA NA 16.1 11.4 2.9 -3.9 NA NA NA NA ...
##  $ total_precip: num  15.5 182.6 57.9 41.7 28.7 ...
##  $ station     : Factor w/ 9 levels "Arkona","Centralia A",..: 1 1 1 1 1 1 8 8 8 8 ...
london_precip <- read.csv("data/london_precip.csv")
london_temp <- read.csv("data/london_temp.csv")
london_precip$month <- as.factor(london_precip$month) #Make a factor 
london_temp$month <- as.factor(london_temp$month) #Make a factor 
climateNA <- read_csv("data/climateNA_monthly.csv") 
## Warning: Missing column names filled in: 'X1' [1]
## Parsed with column specification:
## cols(
##   X1 = col_double(),
##   Year = col_double(),
##   Month = col_double(),
##   Tmax = col_double(),
##   Tmin = col_double(),
##   Tave = col_double(),
##   PPT = col_double(),
##   Rad = col_double(),
##   Eref = col_double(),
##   CMD = col_double(),
##   RH = col_double()
## )
climateNA.temp <- climateNA %>% select(Year, Month, Tave)
climateNA.precip <- climateNA %>% select(Year, Month, PPT)

3 Create Regional Data

3.1 Clean Data

Here I go through each station for the temp and precip data, removing years where all 12 months aren’t present so that I can later calculate the mean values for each year. I did this manually a while ago and there is probably a more efficient way to do it, however it works and it’s already done so I left it as is.

#Arkona
Arkona <- pinery %>% filter(station == "Arkona") #Filter out station of interest

Arkona_T <- Arkona %>% select(-total_precip) %>% filter(!is.na(Arkona$mean_temp)==T) #remove precipitation data and "NA"
table(Arkona_T$year) #Generate table after cleaning
## 
## 1882 1883 1884 1885 1886 1887 1888 1889 1890 1891 1892 1893 1894 1895 1896 1897 
##    4   11   11   12   12   12   12   12   12   12   12    5   12   12   12   12 
## 1898 1899 1900 1901 1902 1903 1904 1905 1906 1907 1908 1909 1910 1911 1912 1913 
##   12   12   10   12    4   12   12   12   12   12   12   12   12   12   12   12 
## 1914 1915 
##   12    3
Arkona_T2 <- Arkona_T %>% filter(year != 1882 & year != 1883 & 
                                   year != 1884 & year != 1893 & year != 1900 & year != 1902 & year != 1915) #remove years with < 12 months
table(Arkona_T2$year) #Generate table after cleaning
## 
## 1885 1886 1887 1888 1889 1890 1891 1892 1894 1895 1896 1897 1898 1899 1901 1903 
##   12   12   12   12   12   12   12   12   12   12   12   12   12   12   12   12 
## 1904 1905 1906 1907 1908 1909 1910 1911 1912 1913 1914 
##   12   12   12   12   12   12   12   12   12   12   12
Arkona_P <- Arkona %>% select(-mean_temp) %>% filter(!is.na(Arkona$total_precip)==T) #remove temperature data and "NA"
#table(Arkona_P$year) #show how many observations in each year
Arkona_P2 <- Arkona_P %>% filter(year != 1882 & year != 1893 & year != 1902 & year != 1915)
#table(Arkona_P2$year) #Generate table after cleaning

#Centralia A
CentraliaA <- pinery %>% filter(station == "Centralia A")

CentraliaA_T <- CentraliaA %>% select(-total_precip) %>% filter(!is.na(CentraliaA$mean_temp)==T) #remove "NA"
#table(CentraliaA_T$year)  
CentraliaA_T2 <- CentraliaA_T %>% filter(year != 1945 & year != 1946 & year != 1947 & year != 1966)
#table(CentraliaA_T2$year)  

CentraliaA_P <- CentraliaA %>% select(-mean_temp) %>% filter(!is.na(CentraliaA$total_precip)==T) #remove "NA"
#table(CentraliaA_P$year)  
CentraliaA_P2 <- CentraliaA_P %>% filter(year != 1946 & year != 1947 & year != 1966)
#table(CentraliaA_P2$year)  

#Exeter
Exeter <- pinery %>% filter(station == "Exeter")

Exeter_T <- Exeter %>% select(-total_precip) %>% filter(!is.na(Exeter$mean_temp)==T) #remove "NA"
#table(Exeter_T$year)  
Exeter_T2 <- Exeter_T %>% filter(year != 1967 & year != 1976 & year != 1985)
#table(Exeter_T2$year)  

Exeter_P <- Exeter %>% select(-mean_temp) %>% filter(!is.na(Exeter$total_precip)==T) #remove "NA"
#table(Exeter_P$year)  
Exeter_P2 <- Exeter_P %>% filter(year != 1961 & year != 1962 & year != 1964 & year != 1976 & year != 1985)
#table(Exeter_P2$year) 

#Forest
Forest <- pinery %>% filter(station == "Forest")

Forest_T <- Forest %>% select(-total_precip) %>% filter(!is.na(Forest$mean_temp)==T) #remove "NA"
#table(Forest_T$year)  
Forest_T2 <- Forest_T %>% filter(year != 1924 & year != 1940 & year != 1962 & year != 1963)
#table(Forest_T2$year) 

Forest_P <- Forest %>% select(-mean_temp) %>% filter(!is.na(Forest$total_precip)==T) #remove "NA"
#table(Forest_P$year)  
Forest_P2 <- Forest_P %>% filter(year > 1933 & year != 1937 & year!= 1943 & year != 1955 & year != 1960 & year < 1960)
#table(Forest_P2$year) 

#Lucan
Lucan <- pinery %>% filter(station == "Lucan")

Lucan_T <- Lucan %>% select(-total_precip) %>% filter(!is.na(Lucan$mean_temp)==T) #remove "NA"
#table(Lucan_T$year)  
Lucan_T2 <- Lucan_T %>% filter(year != 1915 & year != 1924 & year != 1948 & year != 1949 & year != 1950
                               & year != 1951 & year != 1952 & year != 1953 & year != 1954 & year != 1955
                               & year != 1956 & year != 1957 & year != 1958 & year != 1959 & year != 1961)
#table(Lucan_T2$year) 

Lucan_P <- Lucan %>% select(-mean_temp) %>% filter(!is.na(Lucan$total_precip)==T) #remove "NA"
#table(Lucan_P$year)  
Lucan_P2 <- Lucan_P %>% filter(year != 1915 & year != 1923 & year != 1943 & year < 1948)
#table(Lucan_P2$year)

#Pinery
Pinery <- pinery %>% filter(station == "Pinery")

Pinery_T <- Pinery %>% select(-total_precip) %>% filter(!is.na(Pinery$mean_temp)==T) #remove "NA"
#table(Pinery_T$year)  
Pinery_T2 <- Pinery_T %>% filter(year != 1979 & year != 1984)
#table(Pinery_T2$year) 

Pinery_P <- Pinery %>% select(-mean_temp) %>% filter(!is.na(Pinery$total_precip)==T) #remove "NA"
#table(Pinery_P$year)  
Pinery_P2 <- Pinery_P %>% filter(year != 1979 & year != 1984)
#table(Pinery_P2$year)

#Strathrow-Mullifary
Strathroy <- pinery %>% filter(station == "Strathroy-Mullifary")

Strathroy_T <- Strathroy %>% select(-total_precip) %>% filter(!is.na(Strathroy$mean_temp)==T) #remove "NA"
#table(Strathroy_T$year)  
Strathroy_T2 <- Strathroy_T %>% filter(year != 2009 & year != 2010 & year != 2013 & year != 2014 & year != 2015 & year < 2017)
#table(Strathroy_T2$year) 

Strathroy_P <- Strathroy %>% select(-mean_temp) %>% filter(!is.na(Strathroy$total_precip)==T) #remove "NA"
#table(Strathroy_P$year)  
Strathroy_P2 <- Strathroy_P %>% filter(year != 2009 & year != 2010 & year != 2012 & year != 2013 
                                       & year != 2014 & year != 2015 & year < 2017)
#table(Strathroy_P2$year)

#Thedford (no temp)
Thedford <- pinery %>% filter(station == "Thedford")

Thedford_P <- Thedford %>% select(-mean_temp) %>% filter(!is.na(Thedford$total_precip)==T) #remove "NA"
#table(Thedford_P$year)  
Thedford_P2 <- Thedford_P %>% filter(year != 1883 & year != 1897)
#table(Thedford_P2$year)

#Watford (can't use precipitation data)
Watford <- pinery %>% filter(station == "Watford")

Watford_T <- Watford %>% select(-total_precip) %>% filter(!is.na(Watford$mean_temp)==T) #remove "NA"
#table(Watford_T$year)  
Watford_T2 <- Watford_T %>% filter(year == 1927 | year == 1928)
#table(Watford_T2$year) 

Watford_P <- Watford %>% select(-mean_temp) %>% filter(!is.na(Watford$total_precip)==T) #remove "NA"
#table(Watford_P$year)  

Re-combine Data

pinery_temp <- rbind(Arkona_T2, CentraliaA_T2, Exeter_T2, Forest_T2, Lucan_T2, Pinery_T2, Strathroy_T2, Watford_T2)
pinery_precip <- rbind(Arkona_P2, CentraliaA_P2, Exeter_P2, Forest_P2, Lucan_P2, Pinery_P2, Strathroy_P2, Thedford_P2)

Ensure the cleaning worked by plotting histograms of the number of observations (months) in each year (should only be multiples of 12)

#histogram (temperature)
graph <- ggplot(pinery_temp, aes(x = year)) +
  geom_histogram(position = 'identity', binwidth = 1, color = "black", fill = "grey") +
  geom_hline(yintercept = 12, linetype="solid", color = "red") + #shows 12 months
    geom_hline(yintercept = 24, linetype="solid", color = "red") + 
    geom_hline(yintercept = 36, linetype="solid", color = "red") + 
  theme_bw()
graph

#histogram (precipitation)
graph <- ggplot(pinery_precip, aes(x = year)) +
  geom_histogram(position = 'identity', binwidth = 1, color = "black", fill = "grey") +
   geom_hline(yintercept = 12, linetype="solid", color = "red") + #shows 12 months
    geom_hline(yintercept = 24, linetype="solid", color = "red") + 
  theme_bw()
graph

3.2 Create Composite Record

Calculate one mean temperature and total precipitation for each year/month where more than one station of data is available

#Temperature
str(pinery_temp)
## 'data.frame':    1944 obs. of  4 variables:
##  $ year     : int  1885 1885 1885 1885 1885 1885 1885 1885 1885 1885 ...
##  $ month    : int  1 2 3 4 5 6 7 8 9 10 ...
##  $ mean_temp: num  -9.4 -13.5 -8.9 4.2 12.1 16.4 21.1 17 14.8 8 ...
##  $ station  : Factor w/ 9 levels "Arkona","Centralia A",..: 1 1 1 1 1 1 1 1 1 1 ...
pinery_temp_mean <- pinery_temp %>%
  group_by(year,month) %>%
  summarise(mean_temp_year = mean(mean_temp), sd_temp = sd(mean_temp), n = n()) %>%
  mutate(se_temp = sd_temp / sqrt(n), lower = (mean_temp_year - se_temp), upper = (mean_temp_year + se_temp)) %>% round(2)
## `summarise()` regrouping output by 'year' (override with `.groups` argument)
pinery_temp_mean$month <- as.factor(pinery_temp_mean$month)

#Precipitation
pinery_precip_mean <- pinery_precip %>%
  group_by(year,month) %>%
  summarise(total_precip_year = mean(total_precip), sd_precip = sd(total_precip), n = n()) %>%
  mutate(se_precip = sd_precip / sqrt(n), lower = total_precip_year - sd_precip, upper = total_precip_year + sd_precip) %>% round(2)
## `summarise()` regrouping output by 'year' (override with `.groups` argument)
head(pinery_precip_mean)
## # A tibble: 6 x 8
## # Groups:   year [1]
##    year month total_precip_year sd_precip     n se_precip lower upper
##   <dbl> <dbl>             <dbl>     <dbl> <dbl>     <dbl> <dbl> <dbl>
## 1  1883     1              30.2        NA     1        NA    NA    NA
## 2  1883     2              94.7        NA     1        NA    NA    NA
## 3  1883     3              70.6        NA     1        NA    NA    NA
## 4  1883     4              33          NA     1        NA    NA    NA
## 5  1883     5             112.         NA     1        NA    NA    NA
## 6  1883     6             158.         NA     1        NA    NA    NA
ggplot(pinery_precip_mean, aes(x = year)) +
  geom_histogram(position = 'identity', binwidth = 1, color = "black", fill = "grey")

ggplot(pinery_temp_mean, aes(x = year)) +
  geom_histogram(position = 'identity', binwidth = 1, color = "black", fill = "grey")

3.3 Fill Empty Years

3.3.1 Temperature

Fill empty months-years with data from London

p.t <- pinery_temp_mean %>% select(year, month, mean_temp_year)

l.t <- london_temp %>% filter(month == "Jan" |
                                month == "Feb" |
                                month == "Mar" |
                                month == "Apr" |
                                month == "May" |
                                month == "Jun" |
                                month == "Jul" |
                                month == "Aug" |
                                month == "Sep" |
                                month == "Oct" |
                                month == "Nov" |
                                month == "Dec") %>% 
  mutate(month = recode(month, "Jan"="1", 
                                     "Feb"="2", 
                                     "Mar"="3",
                                     "Apr"="4",
                                     "May"="5", 
                                     "Jun"="6", 
                                     "Jul"="7",
                                     "Aug"="8",
                                     "Sep"="9",
                                     "Oct"="10",
                                     "Nov"="11",
                                     "Dec"="12"))
PL.temp <- full_join(l.t, p.t) %>%
  rename("pinery" = "mean_temp_year") %>%
  rename("london" = "mean_temp")
## Joining, by = c("year", "month")
head(PL.temp)
##   year month london pinery
## 1 1883     1     NA     NA
## 2 1884     1  -10.1     NA
## 3 1885     1   -9.2   -9.4
## 4 1886     1   -7.2   -7.2
## 5 1887     1   -8.3   -8.4
## 6 1888     1  -10.1   -9.6

Select all rows with regional data - source = Pinery Select all rows with NA - source = London

PL.1 <- PL.temp %>%  #use london value where pinery is absent
  filter(is.na(pinery)) %>%
  mutate(temp = london)
PL.1$source <- "London"

PL.2 <- PL.temp %>% # use pinery value where pinery is present
  filter(!is.na(pinery)) %>%
  mutate(temp = pinery)
PL.2$source <- "Pinery"

PL.temp2 <- rbind(PL.1, PL.2) %>%
  arrange(year)

(nrow(PL.2) + nrow(PL.1)) == nrow(PL.temp) #check equal observations
## [1] TRUE
PL.temp2 %>% filter(is.na(temp))
##   year month london pinery temp source
## 1 1883     1     NA     NA   NA London
## 2 1883     2     NA     NA   NA London
## 3 2014     9     NA     NA   NA London

Fill leftover NAs with climateNA values

head(climateNA.temp)
## # A tibble: 6 x 3
##    Year Month  Tave
##   <dbl> <dbl> <dbl>
## 1  1901     1  -4.2
## 2  1902     1  -5.6
## 3  1903     1  -5.1
## 4  1904     1  -9.5
## 5  1905     1  -7.5
## 6  1906     1  -0.7
climateNA.temp <- climateNA.temp %>%
  rename("month" = "Month") %>%
  rename("year" = "Year") %>%
  rename("climateNA" = "Tave")

climateNA.temp$month <- as.factor(climateNA.temp$month)

PL.temp3 <- full_join(PL.temp2, climateNA.temp)
## Joining, by = c("year", "month")
PL.temp3 <- PL.temp3[, c(1:4,7,5,6)] #re-order columns

head(PL.temp3)
##   year month london pinery climateNA temp source
## 1 1883     1     NA     NA        NA   NA London
## 2 1883     2     NA     NA        NA   NA London
## 3 1883     3   -6.5     NA        NA -6.5 London
## 4 1883     4    5.1     NA        NA  5.1 London
## 5 1883     5   10.0     NA        NA 10.0 London
## 6 1883     6   17.3     NA        NA 17.3 London
PL.3 <- PL.temp3 %>%  
  filter(is.na(temp)) %>%
  mutate(temp = climateNA)
PL.3$source <- "ClimateNA"

PL.4 <- PL.temp3 %>% 
  filter(!is.na(temp))

PL.temp4 <- rbind(PL.3, PL.4) %>%
  arrange(year) %>% 
  filter(year > 1890) #missing values from 1890s remain

(nrow(PL.3) + nrow(PL.4)) == nrow(PL.temp3) #check equal observations
## [1] TRUE
PL.temp4 %>% filter(is.na(temp)) #values still missing
## [1] year      month     london    pinery    climateNA temp      source   
## <0 rows> (or 0-length row.names)
ggplot(PL.temp4, aes(x = year)) +
  geom_histogram(position = 'identity', binwidth = 1, color = "black", fill = "grey") +
  geom_hline(yintercept = 12, linetype="solid", color = "red") + #shows 12 months
  theme_bw()

3.3.2 Precipitation

p.p <- pinery_precip_mean %>% select(year, month, total_precip_year)

l.p <- london_precip %>% filter(month == "Jan" |
                                month == "Feb" |
                                month == "Mar" |
                                month == "Apr" |
                                month == "May" |
                                month == "Jun" |
                                month == "Jul" |
                                month == "Aug" |
                                month == "Sep" |
                                month == "Oct" |
                                month == "Nov" |
                                month == "Dec") %>% 
  mutate(month = recode(month, "Jan"="1", 
                                     "Feb"="2", 
                                     "Mar"="3",
                                     "Apr"="4",
                                     "May"="5", 
                                     "Jun"="6", 
                                     "Jul"="7",
                                     "Aug"="8",
                                     "Sep"="9",
                                     "Oct"="10",
                                     "Nov"="11",
                                     "Dec"="12"))
p.p$month <- as.factor(p.p$month)
PL.precip <- full_join(l.p, p.p) %>%
  rename("pinery" = "total_precip_year") %>%
  rename("london" = "total_precip")
## Joining, by = c("year", "month")
PL.5 <- PL.precip %>%  #use london value where pinery is absent
  filter(is.na(pinery)) %>%
  mutate(precip = london)
PL.5$source <- "London"

PL.6 <- PL.precip %>% # use pinery value where pinery is present
  filter(!is.na(pinery)) %>% 
  mutate(precip = pinery)
PL.6$source <- "Pinery"

PL.precip2 <- rbind(PL.5, PL.6) %>%
  arrange(year)

(nrow(PL.5) + nrow(PL.6)) == nrow(PL.precip2) #check equal observations
## [1] TRUE
PL.precip2 %>% filter(is.na(precip)) #values still missing
##    year month london pinery precip source
## 1  2009     4     NA     NA     NA London
## 2  2009     5     NA     NA     NA London
## 3  2009     6     NA     NA     NA London
## 4  2009     7     NA     NA     NA London
## 5  2009     8     NA     NA     NA London
## 6  2009     9     NA     NA     NA London
## 7  2009    10     NA     NA     NA London
## 8  2010     4     NA     NA     NA London
## 9  2010     5     NA     NA     NA London
## 10 2010     6     NA     NA     NA London
## 11 2010     7     NA     NA     NA London
## 12 2010     8     NA     NA     NA London
## 13 2010     9     NA     NA     NA London
## 14 2010    10     NA     NA     NA London
## 15 2012     4     NA     NA     NA London
## 16 2012     5     NA     NA     NA London
## 17 2012     6     NA     NA     NA London
## 18 2012     7     NA     NA     NA London
## 19 2012     8     NA     NA     NA London
## 20 2012     9     NA     NA     NA London
## 21 2012    10     NA     NA     NA London
## 22 2013     4     NA     NA     NA London
## 23 2013     5     NA     NA     NA London
## 24 2013     6     NA     NA     NA London
## 25 2013     7     NA     NA     NA London
## 26 2013     8     NA     NA     NA London
## 27 2013     9     NA     NA     NA London
## 28 2013    10     NA     NA     NA London
## 29 2014     4     NA     NA     NA London
## 30 2014     5     NA     NA     NA London
## 31 2014     6     NA     NA     NA London
## 32 2014     7     NA     NA     NA London
## 33 2014     8     NA     NA     NA London
## 34 2014     9     NA     NA     NA London
## 35 2014    10     NA     NA     NA London
## 36 2015     4     NA     NA     NA London
## 37 2015     5     NA     NA     NA London
## 38 2015     6     NA     NA     NA London
## 39 2015     7     NA     NA     NA London
## 40 2015     8     NA     NA     NA London
## 41 2015     9     NA     NA     NA London
## 42 2015    10     NA     NA     NA London
## 43 2017     4     NA     NA     NA London
## 44 2017     5     NA     NA     NA London
## 45 2017     6     NA     NA     NA London
## 46 2017     7     NA     NA     NA London
## 47 2017     8     NA     NA     NA London
## 48 2017     9     NA     NA     NA London
## 49 2017    10     NA     NA     NA London
## 50 2017    11     NA     NA     NA London
## 51 2017    12     NA     NA     NA London

Fill leftover NAs with climateNA values

head(climateNA.temp)
## # A tibble: 6 x 3
##    year month climateNA
##   <dbl> <fct>     <dbl>
## 1  1901 1          -4.2
## 2  1902 1          -5.6
## 3  1903 1          -5.1
## 4  1904 1          -9.5
## 5  1905 1          -7.5
## 6  1906 1          -0.7
climateNA.precip <- climateNA.precip %>%
  rename("month" = "Month") %>%
  rename("year" = "Year") %>%
  rename("climateNA" = "PPT")

climateNA.precip$month <- as.factor(climateNA.precip$month)

PL.precip3 <- full_join(PL.precip2, climateNA.precip)
## Joining, by = c("year", "month")
PL.precip3 <- PL.precip3[, c(1:4,7,5,6)] #re-order columns

head(PL.temp3)
##   year month london pinery climateNA temp source
## 1 1883     1     NA     NA        NA   NA London
## 2 1883     2     NA     NA        NA   NA London
## 3 1883     3   -6.5     NA        NA -6.5 London
## 4 1883     4    5.1     NA        NA  5.1 London
## 5 1883     5   10.0     NA        NA 10.0 London
## 6 1883     6   17.3     NA        NA 17.3 London
PL.7 <- PL.precip3 %>%  
  filter(is.na(precip)) %>%
  mutate(precip = climateNA)
PL.7$source <- "ClimateNA"

PL.8 <- PL.precip3 %>% 
  filter(!is.na(precip))

PL.7
##    year month london pinery climateNA precip    source
## 1  2009     4     NA     NA       132    132 ClimateNA
## 2  2009     5     NA     NA        67     67 ClimateNA
## 3  2009     6     NA     NA       108    108 ClimateNA
## 4  2009     7     NA     NA        72     72 ClimateNA
## 5  2009     8     NA     NA        90     90 ClimateNA
## 6  2009     9     NA     NA        40     40 ClimateNA
## 7  2009    10     NA     NA       116    116 ClimateNA
## 8  2010     4     NA     NA        64     64 ClimateNA
## 9  2010     5     NA     NA       103    103 ClimateNA
## 10 2010     6     NA     NA       132    132 ClimateNA
## 11 2010     7     NA     NA       100    100 ClimateNA
## 12 2010     8     NA     NA        36     36 ClimateNA
## 13 2010     9     NA     NA        76     76 ClimateNA
## 14 2010    10     NA     NA        57     57 ClimateNA
## 15 2012     4     NA     NA        58     58 ClimateNA
## 16 2012     5     NA     NA        50     50 ClimateNA
## 17 2012     6     NA     NA        52     52 ClimateNA
## 18 2012     7     NA     NA        67     67 ClimateNA
## 19 2012     8     NA     NA        70     70 ClimateNA
## 20 2012     9     NA     NA        73     73 ClimateNA
## 21 2012    10     NA     NA       128    128 ClimateNA
## 22 2013     4     NA     NA       115    115 ClimateNA
## 23 2013     5     NA     NA        71     71 ClimateNA
## 24 2013     6     NA     NA       131    131 ClimateNA
## 25 2013     7     NA     NA        91     91 ClimateNA
## 26 2013     8     NA     NA        85     85 ClimateNA
## 27 2013     9     NA     NA        51     51 ClimateNA
## 28 2013    10     NA     NA       123    123 ClimateNA
## 29 2014     4     NA     NA        81     81 ClimateNA
## 30 2014     5     NA     NA        81     81 ClimateNA
## 31 2014     6     NA     NA        88     88 ClimateNA
## 32 2014     7     NA     NA        90     90 ClimateNA
## 33 2014     8     NA     NA       103    103 ClimateNA
## 34 2014     9     NA     NA        95     95 ClimateNA
## 35 2014    10     NA     NA        84     84 ClimateNA
## 36 2015     4     NA     NA        74     74 ClimateNA
## 37 2015     5     NA     NA        95     95 ClimateNA
## 38 2015     6     NA     NA       118    118 ClimateNA
## 39 2015     7     NA     NA        65     65 ClimateNA
## 40 2015     8     NA     NA        79     79 ClimateNA
## 41 2015     9     NA     NA        59     59 ClimateNA
## 42 2015    10     NA     NA        68     68 ClimateNA
## 43 2017     4     NA     NA       118    118 ClimateNA
## 44 2017     5     NA     NA       122    122 ClimateNA
## 45 2017     6     NA     NA       116    116 ClimateNA
## 46 2017     7     NA     NA        75     75 ClimateNA
## 47 2017     8     NA     NA        82     82 ClimateNA
## 48 2017     9     NA     NA        86     86 ClimateNA
## 49 2017    10     NA     NA       107    107 ClimateNA
## 50 2017    11     NA     NA       110    110 ClimateNA
## 51 2017    12     NA     NA        50     50 ClimateNA
## 52 2018     1     NA     NA        66     66 ClimateNA
## 53 2018     2     NA     NA       107    107 ClimateNA
## 54 2018     3     NA     NA        68     68 ClimateNA
## 55 2018     4     NA     NA       102    102 ClimateNA
## 56 2018     5     NA     NA        67     67 ClimateNA
## 57 2018     6     NA     NA        97     97 ClimateNA
## 58 2018     7     NA     NA        89     89 ClimateNA
## 59 2018     8     NA     NA       180    180 ClimateNA
## 60 2018     9     NA     NA        71     71 ClimateNA
## 61 2018    10     NA     NA       100    100 ClimateNA
## 62 2018    11     NA     NA       136    136 ClimateNA
## 63 2018    12     NA     NA        58     58 ClimateNA
PL.precip4 <- rbind(PL.7, PL.8) %>%
  arrange(year) %>%
  filter(year > 1890)

(nrow(PL.7) + nrow(PL.8)) == nrow(PL.precip3) #check equal observations
## [1] TRUE
PL.precip4 %>% filter(is.na(precip)) #values still missing
## [1] year      month     london    pinery    climateNA precip    source   
## <0 rows> (or 0-length row.names)
ggplot(PL.precip4, aes(x = year)) +
  geom_histogram(position = 'identity', binwidth = 1, color = "black", fill = "grey") +
  geom_hline(yintercept = 12, linetype="solid", color = "red") + #shows 12 months
  theme_bw()

3.3.3 Combine

Fix levels

PL.temp4 <- droplevels(PL.temp4) #drops un-used levels
PL.temp4$month <- factor(PL.temp4$month, levels = c("1","2","3","4", "5", "6", "7", "8", "9", "10", "11", "12")) #re-level
levels(PL.temp4$month)
##  [1] "1"  "2"  "3"  "4"  "5"  "6"  "7"  "8"  "9"  "10" "11" "12"
PL.precip4 <- droplevels(PL.precip4) #drops un-used levels
PL.precip4$month <- factor(PL.precip4$month, levels = c("1","2","3","4", "5", "6", "7", "8", "9", "10", "11", "12")) #re-level
levels(PL.precip4$month)
##  [1] "1"  "2"  "3"  "4"  "5"  "6"  "7"  "8"  "9"  "10" "11" "12"
PL.temp.only <- PL.temp4 %>% select(year, month, temp)
PL.precip.only <- PL.precip4 %>% select(year, month, precip)

PL.clean <- full_join(PL.temp.only, PL.precip.only, by = c("year", "month"))

str(PL.clean)
## 'data.frame':    1536 obs. of  4 variables:
##  $ year  : num  1891 1891 1891 1891 1891 ...
##  $ month : Factor w/ 12 levels "1","2","3","4",..: 1 2 3 4 5 6 7 8 9 10 ...
##  $ temp  : num  -4.2 -1.9 -2.5 7 10.9 18.8 17.5 18.8 16.9 8.6 ...
##  $ precip: num  53.6 72.5 89 45.4 40 ...
PL.clean %>% filter(is.na(temp)) #values still missing
## [1] year   month  temp   precip
## <0 rows> (or 0-length row.names)
nrow(PL.clean) == (2019-1891)*12 #check 12 months / year for the 128 years
## [1] TRUE

4 Monthly Data

4.1 Climograph

study_t <- PL.clean %>% 
  group_by(month) %>% 
  summarize(mean_temp = mean(temp) %>% round(2)) 
## `summarise()` ungrouping output (override with `.groups` argument)
head(study_t)
## # A tibble: 6 x 2
##   month mean_temp
##   <fct>     <dbl>
## 1 1         -5.68
## 2 2         -5.68
## 3 3         -0.46
## 4 4          6.38
## 5 5         12.8 
## 6 6         18.2
study_p <- PL.clean %>%
  group_by(month) %>%
  summarize(total_precip= mean(precip), sd_precip = sd(precip), n = n()) %>%
  mutate(se_precip = sd_precip / sqrt(n), lower = total_precip - se_precip, upper = total_precip + se_precip %>% round(2)) 
## `summarise()` ungrouping output (override with `.groups` argument)
head(study_p)
## # A tibble: 6 x 7
##   month total_precip sd_precip     n se_precip lower upper
##   <fct>        <dbl>     <dbl> <int>     <dbl> <dbl> <dbl>
## 1 1             82.8      33.3   128      2.94  79.8  85.7
## 2 2             67.3      28.5   128      2.52  64.8  69.8
## 3 3             65.7      28.3   128      2.50  63.2  68.2
## 4 4             74.6      32.0   128      2.83  71.8  77.5
## 5 5             81.5      36.9   128      3.26  78.2  84.8
## 6 6             81.0      37.7   128      3.33  77.6  84.3
p.climograph <- ggplot() +
  geom_boxplot(data = PL.clean, aes(x = month, y = temp), fill = "red3", alpha = 0.7) +
  geom_path(data = study_p, aes(x = month, y = total_precip-75, group=1), col= "steelblue4", size = 1) +
  geom_ribbon(data = study_p, aes(x = month, ymin = (lower-75), ymax = (upper-75), group=1), fill = "steelblue4", alpha = 0.2) +
  scale_y_continuous("Temperature (°C)", sec.axis = sec_axis(~ . + 75, name = "Precipitation (mm)")) +
  scale_x_discrete("Month") +
  theme_bw() +
   theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())

p.climograph

#ggsave('p.climograph.pdf', p.climograph, units = 'cm', width = 20, height = 12)
ggplot(data = PL.temp4, aes(x=year, y = temp)) +
  geom_point(aes(col=source), size=0.3) +
  geom_line(alpha = 0.6) +
  facet_grid(month~.) + 
  ylab("Mean Monthly Temperature (°C)") +
  theme_classic() 

ggplot(data = PL.temp4, aes(x=year, y = temp)) +
  geom_point(aes(col=source), size=0.3) +
  facet_grid(month~source) + 
  ylab("Mean Monthly Temperature (°C)") +
  theme_classic() 

ggplot(data = PL.precip4, aes(x=year, y = precip)) +
  geom_point(aes(col=source), size=0.3) +
  geom_line(alpha = 0.6) +
  facet_grid(month~.) + 
  ylab("Total Monthly Precipitation (mm)") +
  theme_classic() 

ggplot(data = PL.precip4, aes(x=year, y = precip)) +
  geom_point(aes(col=source), size=0.3) +
  facet_grid(month~source) + 
  ylab("Total Monthly Precipitation (mm)") +
  theme_classic() 

ggplot(data = PL.temp4, aes(x = month, y = temp)) + 
  geom_boxplot(alpha = 0.7, outlier.shape=NA) + 
  geom_jitter(aes(col=source), size=0.6, alpha=0.5) +
  ylab("Mean Monthly Temperature (°C)") +
  theme_bw() +
   theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) 

ggplot(data = PL.precip4, aes(x = month, y = precip)) + 
  geom_boxplot(alpha = 0.7, outlier.shape=NA) + 
  geom_jitter(aes(col=source), size=0.6, alpha=0.5) +
  ylab("Total Monthly Precipitation (mm)") +
  theme_bw() +
   theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) 

5 Annual Data

5.1 Clean Record

clean_temp_annual <- PL.clean %>% group_by(year) %>% summarise(mean_temp = mean(temp))
## `summarise()` ungrouping output (override with `.groups` argument)
clean_precip_annual <- PL.clean %>% group_by(year) %>% summarise(mean_precip = mean(precip))
## `summarise()` ungrouping output (override with `.groups` argument)
ggplot() +
  geom_line(data=clean_temp_annual, aes(x = year, y = mean_temp), colour="lightsalmon4") +
  #Aesthetics
  xlab("Year") +
  ylab("Mean Annual Temperature (°C)") +
  scale_x_continuous(limits = c(1890, 2020), breaks=seq(1900,2020,20)) +
  theme_bw() + theme(plot.title = element_text(hjust = 0.5)) +
 theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) 

ggplot() +
  geom_line(data=clean_precip_annual, aes(x = year, y = mean_precip), colour="lightsalmon4") +
  xlab("Year") +
  ylab("Total Annual Precipitation (mm)") +
  scale_x_continuous(limits = c(1890, 2020), breaks=seq(1900,2020,20)) +
 theme_bw() + theme(plot.title = element_text(hjust = 0.5)) +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) 

5.2 Individual Records

london_temp_annual <- london_temp %>% filter(month == "Annual")
london_precip_annual <- london_precip %>% filter(month == "Annual")
pinery_temp_annual <- pinery_temp_mean %>% group_by(year) %>% summarise(mean_temp = mean(mean_temp_year))
## `summarise()` ungrouping output (override with `.groups` argument)
pinery_precip_annual <- pinery_precip_mean %>% group_by(year) %>% summarise(total_precip = sum(total_precip_year))
## `summarise()` ungrouping output (override with `.groups` argument)
london_temp_2 <- london_temp_annual %>% select(-month)
london_temp_2$location <- "London"
pinery_temp_annual$location <- "Pinery"
annual_temp <- rbind(london_temp_2,pinery_temp_annual) #join data sets

london_precip_2 <- london_precip_annual %>% select(-month)
london_precip_2$location <- "London"
pinery_precip_annual$location <- "Pinery"
annual_precip <- rbind(london_precip_2,pinery_precip_annual) #join data sets
annual_temp_plot <- ggplot() +
  # London Data
  geom_line(data=london_temp_annual, aes(x = year, y = mean_temp), colour = "lightseagreen", linetype = "dashed", alpha = 0.8) +
  # Pinery Data
  geom_line(data=pinery_temp_annual, aes(x = year, y = mean_temp), colour="lightsalmon4") +
  #Aesthetics
  xlab("Year") +
  ylab("Mean Annual Temperature (°C)") +
  scale_x_continuous(limits = c(1880, 2020), breaks=seq(1880,2020,20)) +
  #ggtitle("Average Annual Temperature in \n Southwestern Ontario from 1883 to 2017") + 
  theme_bw() + theme(plot.title = element_text(hjust = 0.5)) +
 theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) 

annual_temp_plot
## Warning: Removed 1 row(s) containing missing values (geom_path).

#ggsave('annual_temp.pdf', annual_temp_plot, units = 'cm', width = 15, height = 10)

annual_precip_plot <- ggplot() +
  # London Data
  geom_line(data=london_precip_annual, aes(x = year, y = total_precip), colour = "lightseagreen", linetype = "dashed", alpha = 0.8) +
   # geom_smooth(data=london_precip_annual, aes(x = year, y = total_precip), se = TRUE, 
             # method = 'lm', size = 0.8, colour = "black", fill = "blue") +
  # Pinery Data
  geom_line(data=pinery_precip_annual, aes(x = year, y = total_precip, group = 1),  colour="lightsalmon4")+
   # geom_smooth(data=pinery_precip_annual, aes(x = year, y = total_precip), se = TRUE, 
            #  method = 'lm', size = 0.8, colour = "black", fill = "darkgreen") +
  #Aesthetics
  xlab("Year") +
  ylab("Total Annual Precipitation (mm)") +
  scale_x_continuous(breaks=seq(1880,2020,20)) +
  #ggtitle("Annual Precipitation in Southwestern Ontario \nfrom 1883 to 2017") + 
  scale_colour_manual(name="Location", breaks = c("London", "Pinery"),
                      values = c("blue","darkgreen")) +
  theme_bw() + theme(plot.title = element_text(hjust = 0.5)) +
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank()) 

annual_precip_plot
## Warning: Removed 17 row(s) containing missing values (geom_path).

#ggsave('annual_precip.pdf', annual_precip_plot, units = 'cm', width = 15, height = 10)

annual_plot <- grid.arrange(annual_temp_plot, annual_precip_plot, ncol= 2)
## Warning: Removed 1 row(s) containing missing values (geom_path).
## Warning: Removed 17 row(s) containing missing values (geom_path).

#ggsave('annual_plot.pdf', annual_plot, units = 'cm', width = 23, height = 8)

6 Seasonal Data

winter_t <- pinery_temp %>% filter(month == 1 | month == 2 | month == 3)
winter_t$season <- "winter"
spring_t <- pinery_temp %>% filter(month == 4 | month == 5 | month == 6)
spring_t$season <- "spring"
summer_t <- pinery_temp %>% filter(month == 7 | month == 8 | month == 9)
summer_t$season <- "summer"
fall_t <- pinery_temp %>% filter(month == 10 | month == 11 | month == 12)
fall_t$season <- "fall"
pinery_temp <- rbind(winter_t, spring_t, summer_t, fall_t)
pinery_temp$season <- as.factor(pinery_temp$season) #Make a factor 
str(pinery_temp)
## 'data.frame':    1944 obs. of  5 variables:
##  $ year     : int  1885 1885 1885 1886 1886 1886 1887 1887 1887 1888 ...
##  $ month    : int  1 2 3 1 2 3 1 2 3 1 ...
##  $ mean_temp: num  -9.4 -13.5 -8.9 -7.2 -6.9 -1 -8.4 -5.4 -3.8 -9.6 ...
##  $ station  : Factor w/ 9 levels "Arkona","Centralia A",..: 1 1 1 1 1 1 1 1 1 1 ...
##  $ season   : Factor w/ 4 levels "fall","spring",..: 4 4 4 4 4 4 4 4 4 4 ...
pinery_temp_station <- pinery_temp %>% 
 group_by(year, season, station) %>% 
 summarise(mean_temp = mean(mean_temp))
## `summarise()` regrouping output by 'year', 'season' (override with `.groups` argument)
pinery_temp_season <- pinery_temp_station %>%
  group_by(year, season) %>%
  summarise(mean_temp_season = mean(mean_temp) %>% round(2), sd_temp = sd(mean_temp) %>% round(2), n = n()) %>% 
  mutate(se_temp = (sd_temp / sqrt(n)) %>% round(2), lower = mean_temp_season - se_temp, upper = mean_temp_season + se_temp) 
## `summarise()` regrouping output by 'year' (override with `.groups` argument)
winter_p <- pinery_precip %>% filter(month == 1 | month == 2 | month == 3)
winter_p$season <- "winter"
spring_p <- pinery_precip %>% filter(month == 4 | month == 5 | month == 6)
spring_p$season <- "spring"
summer_p <- pinery_precip %>% filter(month == 7 | month == 8 | month == 9)
summer_p$season <- "summer"
fall_p <- pinery_precip %>% filter(month == 10 | month == 11 | month == 12)
fall_p$season <- "fall"
pinery_precip <- rbind(winter_p, spring_p, summer_p, fall_p)
pinery_precip$season <- as.factor(pinery_precip$season) #Make a factor 

pinery_precip_station <- pinery_precip %>%
  group_by(year, season, station) %>%
  summarise(total_precip = mean(total_precip))
## `summarise()` regrouping output by 'year', 'season' (override with `.groups` argument)
pinery_precip_season <- pinery_precip_station %>% 
  group_by(year, season) %>% 
  summarise(total_precip_season = mean(total_precip) %>% round(2), sd_precip = sd(total_precip) %>% round(2), n = n()) %>%
  mutate(se_precip = sd_precip / sqrt(n), lower = total_precip_season - sd_precip, upper = total_precip_season + sd_precip %>% round(2)) 
## `summarise()` regrouping output by 'year' (override with `.groups` argument)
winter_tl <- london_temp %>% filter(month == "Jan" | month == "Feb" | month == "Mar")
winter_tl$season <- "winter"
spring_tl <- london_temp %>% filter(month == "Apr" | month == "May" | month == "Jun")
spring_tl$season <- "spring"
summer_tl <- london_temp %>% filter(month == "Jul" | month == "Aug" | month == "Sep")
summer_tl$season <- "summer"
fall_tl <- london_temp %>% filter(month == "Oct" | month == "Nov" | month == "Dec")
fall_tl$season <- "fall"
london_temp <- rbind(winter_tl, spring_tl, summer_tl, fall_tl)
london_temp$season <- as.factor(london_temp$season) #Make a factor 

london_temp_season <- london_temp %>% 
  group_by(year, season) %>% 
  summarise(mean_temp_season = mean(mean_temp))
## `summarise()` regrouping output by 'year' (override with `.groups` argument)
winter_pl <- london_precip %>% filter(month == "Jan" | month == "Feb" | month == "Mar")
winter_pl$season <- "winter"
spring_pl <- london_precip %>% filter(month == "Apr" | month == "May" | month == "Jun")
spring_pl$season <- "spring"
summer_pl <- london_precip %>% filter(month == "Jul" | month == "Aug" | month == "Sep")
summer_pl$season <- "summer"
fall_pl <- london_precip %>% filter(month == "Oct" | month == "Nov" | month == "Dec")
fall_pl$season <- "fall"
london_precip <- rbind(winter_pl, spring_pl, summer_pl, fall_pl)
london_precip$season <- as.factor(london_precip$season) #Make a factor 

london_precip_season <- london_precip %>% 
  group_by(year, season) %>% 
  summarise(total_precip_season = mean(total_precip))
## `summarise()` regrouping output by 'year' (override with `.groups` argument)
london_temp_season$season <- factor(london_temp_season$season, levels = c("winter", "spring", "summer", "fall"))
pinery_temp_season$season <- factor(pinery_temp_season$season, levels = c("winter", "spring", "summer", "fall"))
pinery_temp_season$location <- "Pinery"
london_temp_season$location <- "London"
season_temp <- rbind(pinery_temp_season, london_temp_season)

pinery_precip_season$location <- "Pinery"
london_precip_season$location <- "London"
season_precip <- rbind(pinery_precip_season, london_precip_season)

season_temp$location <- factor(season_temp$location, levels = c("Pinery", "London"))
season_precip$location <- factor(season_precip$location, levels = c("Pinery", "London"))

6.1 Temperature

w_temp <- ggplot() +
  geom_line(data=(season_temp %>% filter(season == "winter")), aes(x = year, y = mean_temp_season, col = location, linetype = location)) +
  #Aesthetics
  xlab("Year") + ylab("Temperature (°C)") + ggtitle("Winter") +
  scale_x_continuous(limits = c(1880, 2020), breaks=seq(1880,2020,20)) +
  scale_color_manual(values = c("lightsalmon4", "lightseagreen")) +
  theme_bw() + theme(plot.title = element_text(hjust = 0.5)) +
 theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), legend.position = "none") 

sp_temp <- ggplot() +
  geom_line(data=(season_temp %>% filter(season == "spring")), aes(x = year, y = mean_temp_season, col = location, linetype = location)) +
  #Aesthetics
  xlab("Year") + ylab("Temperature (°C)") + ggtitle("Spring") +
  scale_x_continuous(limits = c(1880, 2020), breaks=seq(1880,2020,20)) +
    scale_color_manual(values = c("lightsalmon4", "lightseagreen")) +
  theme_bw() + theme(plot.title = element_text(hjust = 0.5)) +
 theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), legend.position = "none") 

su_temp <- ggplot() +
  geom_line(data=(season_temp %>% filter(season == "summer")), aes(x = year, y = mean_temp_season, col = location, linetype = location)) +
  #Aesthetics
  xlab("Year") + ylab("Temperature (°C)") + ggtitle("Summer") +
  scale_x_continuous(limits = c(1880, 2020), breaks=seq(1880,2020,20)) +
    scale_color_manual(values = c("lightsalmon4", "lightseagreen")) +
  theme_bw() + theme(plot.title = element_text(hjust = 0.5)) +
 theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), legend.position = "none") 

fa_temp <- ggplot() +
  geom_line(data=(season_temp %>% filter(season == "fall")), aes(x = year, y = mean_temp_season, col = location, linetype = location)) +
  #Aesthetics
  xlab("Year") + ylab("Temperature (°C)") + ggtitle("Fall") +
  scale_x_continuous(limits = c(1880, 2020), breaks=seq(1880,2020,20)) +
    scale_color_manual(values = c("lightsalmon4", "lightseagreen")) +
  theme_bw() + theme(plot.title = element_text(hjust = 0.5)) +
 theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),  legend.position = "none") 

seasonal_temp_plot <- grid.arrange(w_temp, sp_temp, su_temp, fa_temp, nrow = 2, ncol = 2)
## Warning: Removed 1 row(s) containing missing values (geom_path).

#ggsave('seasonal_temperature.pdf', seasonal_temp_plot, units = 'cm', width = 20, height = 12)

6.2 Precipitation

w_precip <- ggplot() +
  geom_line(data=(season_precip %>% filter(season == "winter")), aes(x = year, y = total_precip_season, col = location, linetype = location)) +
  #Aesthetics
  xlab("Year") + ylab("Precipitation (mm)") + ggtitle("Winter") +
  scale_x_continuous(limits = c(1880, 2020), breaks=seq(1880,2020,20)) +
      scale_color_manual(values = c("lightsalmon4", "lightseagreen")) +
  theme_bw() + theme(plot.title = element_text(hjust = 0.5)) +
 theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), legend.position = "none") 

sp_precip <- ggplot() +
  geom_line(data=(season_precip %>% filter(season == "spring")), aes(x = year, y = total_precip_season, col = location, linetype = location)) +
  #Aesthetics
  xlab("Year") + ylab("Precipitation (mm)") + ggtitle("Spring") +
  scale_x_continuous(limits = c(1880, 2020), breaks=seq(1880,2020,20)) +
      scale_color_manual(values = c("lightsalmon4", "lightseagreen")) +
  theme_bw() + theme(plot.title = element_text(hjust = 0.5)) +
 theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), legend.position = "none") 

su_precip <- ggplot() +
  geom_line(data=(season_precip %>% filter(season == "summer")), aes(x = year, y = total_precip_season, col = location, linetype = location)) +
  #Aesthetics
  xlab("Year") + ylab("Precipitation (mm)") + ggtitle("Summer") +
  scale_x_continuous(limits = c(1880, 2020), breaks=seq(1880,2020,20)) +
      scale_color_manual(values = c("lightsalmon4", "lightseagreen")) +
  theme_bw() + theme(plot.title = element_text(hjust = 0.5)) +
 theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), legend.position = "none") 

fa_precip <- ggplot() +
  geom_line(data=(season_precip %>% filter(season == "fall")), aes(x = year, y = total_precip_season, col = location, linetype = location)) +
  #Aesthetics
  xlab("Year") + ylab("Precipitation (mm)") + ggtitle("Fall") +
  scale_x_continuous(limits = c(1880, 2020), breaks=seq(1880,2020,20)) +
      scale_color_manual(values = c("lightsalmon4", "lightseagreen")) +
  theme_bw() + theme(plot.title = element_text(hjust = 0.5)) +
 theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),  legend.position = "none") 

seasonal_precip_plot <- grid.arrange(w_precip, sp_precip, su_precip, fa_precip, nrow = 2, ncol = 2)
## Warning: Removed 1 row(s) containing missing values (geom_path).
## Warning: Removed 15 row(s) containing missing values (geom_path).

## Warning: Removed 15 row(s) containing missing values (geom_path).
## Warning: Removed 16 row(s) containing missing values (geom_path).

#ggsave('seasonal_precipitation.pdf', seasonal_precip_plot, units = 'cm', width = 20, height = 12)

6.3 Temperature Anomaly

https://tamino.wordpress.com/2018/07/29/why-use-temperature-anomaly/

Average defined as 1951-1980 (NASA standard)

pinery_average <- season_temp %>% filter(location == "Pinery") %>% filter(year > 1950 & year < 1981)
pinery_average <- pinery_average %>% group_by(season) %>% summarise(pinery_average = mean(mean_temp_season))
## `summarise()` ungrouping output (override with `.groups` argument)
london_average <- season_temp %>% filter(location == "London") %>% filter(year > 1950 & year < 1981)
london_average <- london_average %>% group_by(season) %>% summarise(london_average = mean(mean_temp_season))
## `summarise()` ungrouping output (override with `.groups` argument)
pinery_anomaly <- right_join(season_temp, pinery_average)
## Joining, by = "season"
pinery_anomaly <- mutate(pinery_anomaly, anomaly = mean_temp_season - pinery_average) %>% filter(location == "Pinery")

london_anomaly <- right_join(season_temp, london_average)
## Joining, by = "season"
london_anomaly <- mutate(london_anomaly, anomaly = mean_temp_season - london_average) %>% filter(location == "London")

season_anomaly <- full_join(pinery_anomaly, london_anomaly) %>% select(-pinery_average, -london_average)
## Joining, by = c("year", "season", "mean_temp_season", "sd_temp", "n", "se_temp", "lower", "upper", "location", "anomaly")
w_anomaly <- ggplot() +
  geom_line(data=(season_anomaly %>% filter(season == "winter")), aes(x = year, y = anomaly, col = location, linetype = location)) +
  #Aesthetics
  xlab("Year") + ylab("Temperature Anomaly (°C)") + ggtitle("Winter") +
  scale_x_continuous(limits = c(1880, 2020), breaks=seq(1880,2020,20)) +
  theme_bw() + theme(plot.title = element_text(hjust = 0.5)) +
 theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), legend.position = "none") 

sp_anomaly <- ggplot() +
  geom_line(data=(season_anomaly %>% filter(season == "spring")), aes(x = year, y = anomaly, col = location, linetype = location)) +
  #Aesthetics
  xlab("Year") + ylab("Temperature Anomaly (°C)") + ggtitle("Spring") +
  scale_x_continuous(limits = c(1880, 2020), breaks=seq(1880,2020,20)) +
  theme_bw() + theme(plot.title = element_text(hjust = 0.5)) +
 theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), legend.position = "none") 

su_anomaly <- ggplot() +
  geom_line(data=(season_anomaly %>% filter(season == "summer")), aes(x = year, y = anomaly, col = location, linetype = location)) +
  #Aesthetics
  xlab("Year") + ylab("Temperature Anomaly (°C)") + ggtitle("Summer") +
  scale_x_continuous(limits = c(1880, 2020), breaks=seq(1880,2020,20)) +
  theme_bw() + theme(plot.title = element_text(hjust = 0.5)) +
 theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(), legend.position = "none") 

fa_anomaly <- ggplot() +
  geom_line(data=(season_anomaly %>% filter(season == "fall")), aes(x = year, y = anomaly, col = location, linetype = location)) +
  #Aesthetics
  xlab("Year") + ylab("Temperature Anomaly (°C)") + ggtitle("Fall") +
  scale_x_continuous(limits = c(1880, 2020), breaks=seq(1880,2020,20)) +
  theme_bw() + theme(plot.title = element_text(hjust = 0.5)) +
 theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),  legend.position = "none") 

seasonal_anomaly_plot <- grid.arrange(w_anomaly, sp_anomaly, su_anomaly, fa_anomaly, nrow = 2, ncol = 2)
## Warning: Removed 1 row(s) containing missing values (geom_path).

#ggsave('seasonal_anomaly.pdf', seasonal_anomaly_plot, units = 'cm', width = 20, height = 12)

7 Save Outputs

write.csv(PL.clean, "outputs/full_climate.csv")
write.csv(PL.temp4, "outputs/full_temp.csv")
write.csv(PL.precip4, "outputs/full_precip.csv")
write.csv(pinery_temp_mean, "outputs/pinery_monthly_temperature.csv")
write.csv(pinery_precip_mean, "outputs/pinery_monthly_precipitation.csv")
write.csv(season_temp, file = "outputs/pinery_season_temperature.csv")
write.csv(season_precip, file = "outputs/pinery_season_precipitation.csv")
write.csv(pinery_temp_annual, file = "outputs/pinery_annual_temp.csv")

8 Reproducibility

Sys.time()
## [1] "2020-08-06 16:25:48 PDT"
git2r::repository()
## Local:    master /Users/JenBaron/Documents/UWO/UWO NSERC/Statistical Analysis/Climate/Climate_Pinery
## Remote:   master @ origin (https://github.com/JenBaron/Climate_Pinery.git)
## Head:     [78f8528] 2020-08-05: filled .rmd
sessionInfo()
## R version 4.0.0 (2020-04-24)
## Platform: x86_64-apple-darwin17.0 (64-bit)
## Running under: macOS Catalina 10.15.6
## 
## Matrix products: default
## BLAS:   /Library/Frameworks/R.framework/Versions/4.0/Resources/lib/libRblas.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/4.0/Resources/lib/libRlapack.dylib
## 
## locale:
## [1] en_CA.UTF-8/en_CA.UTF-8/en_CA.UTF-8/C/en_CA.UTF-8/en_CA.UTF-8
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] gridExtra_2.3   GGally_1.5.0    forcats_0.5.0   stringr_1.4.0  
##  [5] dplyr_1.0.0     purrr_0.3.4     readr_1.3.1     tidyr_1.0.2    
##  [9] tibble_3.0.1    tidyverse_1.3.0 ggplot2_3.3.1  
## 
## loaded via a namespace (and not attached):
##  [1] tidyselect_1.1.0   xfun_0.14          haven_2.2.0        lattice_0.20-41   
##  [5] colorspace_1.4-1   vctrs_0.3.1        generics_0.0.2     htmltools_0.4.0   
##  [9] yaml_2.2.1         utf8_1.1.4         rlang_0.4.6        pillar_1.4.4      
## [13] glue_1.4.1         withr_2.2.0        DBI_1.1.0          RColorBrewer_1.1-2
## [17] dbplyr_1.4.3       modelr_0.1.6       readxl_1.3.1       plyr_1.8.6        
## [21] lifecycle_0.2.0    munsell_0.5.0      gtable_0.3.0       cellranger_1.1.0  
## [25] rvest_0.3.5        evaluate_0.14      labeling_0.3       knitr_1.28        
## [29] fansi_0.4.1        broom_0.5.6        Rcpp_1.0.4.6       scales_1.1.1      
## [33] backports_1.1.7    jsonlite_1.6.1     farver_2.0.3       fs_1.4.1          
## [37] hms_0.5.3          digest_0.6.25      stringi_1.4.6      grid_4.0.0        
## [41] cli_2.0.2          tools_4.0.0        magrittr_1.5       crayon_1.3.4      
## [45] pkgconfig_2.0.3    ellipsis_0.3.1     xml2_1.3.2         reprex_0.3.0      
## [49] lubridate_1.7.8    reshape_0.8.8      assertthat_0.2.1   rmarkdown_2.1     
## [53] httr_1.4.1         rstudioapi_0.11    R6_2.4.1           git2r_0.26.1      
## [57] nlme_3.1-147       compiler_4.0.0